{
 "cells": [
  {
   "cell_type": "code",
   "id": "initial_id",
   "metadata": {
    "collapsed": true,
    "jupyter": {
     "is_executing": true
    }
   },
   "source": [
    "from scipy.integrate import odeint\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "\n",
    "plt.rcParams['axes.unicode_minus'] = False  # 解决负号显示问题\n",
    "plt.rc('font', size=10)\n",
    "plt.rc('font', family='SimHei')"
   ],
   "execution_count": null,
   "outputs": []
  },
  {
   "metadata": {
    "jupyter": {
     "is_executing": true
    }
   },
   "cell_type": "code",
   "source": [
    "m_A = 4866\n",
    "m_B = 2433\n",
    "f = 4890\n",
    "delta_m = 1165.992\n",
    "omiga = 2.2143\n",
    "k_3 = 167.8395  # 兴波\n",
    "k_5 = 80000  # 弹簧劲度系数\n",
    "k_6 = 10000  # 阻尼系数\n",
    "ro = 1025\n",
    "g = 9.8\n",
    "s = np.pi"
   ],
   "id": "46ec89d2bacd559c",
   "execution_count": null,
   "outputs": []
  },
  {
   "metadata": {
    "jupyter": {
     "is_executing": true
    }
   },
   "cell_type": "code",
   "source": [
    "# 定义一个方程组（微分方程组）\n",
    "def pfun(y, t):\n",
    "    y0, y1, y2, y3 = y\n",
    "    return np.array([\n",
    "        1 / (m_A + delta_m) * (-k_3 * y0 - ro * g * s * y1 + f * np.cos(omiga * t) + k_5 * (y3 - y1) + k_6 * (y2 - y0)),\n",
    "        y0,\n",
    "        1 / m_B * (-k_5 * (y3 - y1) - k_6 * (y2 - y0) ),\n",
    "        y2\n",
    "    ])"
   ],
   "id": "f790812a4ee5d2af",
   "execution_count": null,
   "outputs": []
  },
  {
   "metadata": {
    "jupyter": {
     "is_executing": true
    }
   },
   "cell_type": "code",
   "source": [
    "t = np.arange(0, 28.56122813280971, 0.2)  # 创建自变量序列\n",
    "soli = odeint(pfun, [0, 0, 0, 0], t)  # 求数值解\n",
    "\n",
    "plt.figure(figsize=(5, 2))\n",
    "plt.plot(t, soli[:, 1], 'b', label=\"浮子位移\")\n",
    "plt.plot(t, soli[:, 3], 'r', label=\"振子位移\")\n",
    "plt.legend()\n",
    "plt.show()\n",
    "\n",
    "print(soli[:, 1][50])\n",
    "# print(soli[:, 1][100])\n",
    "# print(soli[:, 1][200])\n",
    "# print(soli[:, 1][300])"
   ],
   "id": "a521ef69dadb1a4a",
   "execution_count": null,
   "outputs": []
  },
  {
   "metadata": {
    "ExecuteTime": {
     "end_time": "2025-08-27T13:46:40.239314900Z",
     "start_time": "2025-08-27T09:56:15.936910Z"
    }
   },
   "cell_type": "code",
   "source": [
    "# # 定义一个方程组（微分方程组）\n",
    "# def pfun2(y, t):\n",
    "#     y0, y1, y2, y3 = y\n",
    "#     return np.array([\n",
    "#         1 / (m_A + delta_m) * (-k_3 * y0 - ro * g * s * y1 + f * np.cos(omiga * t) + k_5 * (y3 - y1) + 10000*(abs(y2-y0))**(0.5) * (y2 - y0)),\n",
    "#         y0,\n",
    "#         1 / m_B * (-k_5 * (y3 - y1) - 10000*(abs(y2-y0))**(0.5) * (y2 - y0) ),\n",
    "#         y2\n",
    "#     ])"
   ],
   "id": "a17b429530f725c5",
   "execution_count": 9,
   "outputs": []
  },
  {
   "metadata": {
    "ExecuteTime": {
     "end_time": "2025-08-27T13:46:40.240317900Z",
     "start_time": "2025-08-27T09:56:16.745366Z"
    }
   },
   "cell_type": "code",
   "source": [
    "# t = np.arange(0, 228.56122813280971, 0.2)  # 创建自变量序列\n",
    "# soli = odeint(pfun2, [0, 0, 0, 0], t)  # 求数值解\n",
    "#\n",
    "# plt.figure(figsize=(5, 2))\n",
    "# plt.plot(t, soli[:, 1], 'b', label=\"浮子位移\")\n",
    "# plt.plot(t, soli[:, 3], 'r', label=\"振子位移\")\n",
    "# plt.legend()\n",
    "# plt.show()\n",
    "#\n",
    "# print(soli[:, 1][50])\n",
    "# print(soli[:, 1][100])\n",
    "# print(soli[:, 1][200])\n",
    "# print(soli[:, 1][300])"
   ],
   "id": "6af30fd6cb923f79",
   "execution_count": 10,
   "outputs": []
  },
  {
   "metadata": {
    "ExecuteTime": {
     "end_time": "2025-08-27T13:46:40.250278500Z",
     "start_time": "2025-08-27T09:56:18.165481Z"
    }
   },
   "cell_type": "code",
   "source": [
    "# 在SA函数外部初始化记录容器\n",
    "history = {\n",
    "    'temperature': [],\n",
    "    'best_P': [],\n",
    "    'current_P': [],\n",
    "    'k': [],\n",
    "}"
   ],
   "id": "fa504dcd8171dd89",
   "execution_count": 11,
   "outputs": []
  },
  {
   "metadata": {
    "ExecuteTime": {
     "end_time": "2025-08-27T13:46:40.252278100Z",
     "start_time": "2025-08-27T09:57:00.272955Z"
    }
   },
   "cell_type": "code",
   "source": [
    "def plot_SA_history():\n",
    "    plt.figure(figsize=(15, 4))\n",
    "\n",
    "    # 子图1：损失函数变化\n",
    "    plt.subplot(1, 3, 1)\n",
    "    plt.plot(history['best_P'], 'r-', label='Best P')\n",
    "    plt.plot(history['current_P'], 'b--', alpha=0.5, label='Current P')\n",
    "    plt.xlabel('Iteration')\n",
    "    plt.ylabel('P')\n",
    "    plt.title('P变化曲线')\n",
    "    plt.legend()\n",
    "    plt.grid(True, which=\"both\", ls=\"--\")\n",
    "\n",
    "    # 子图2：温度衰减曲线\n",
    "    plt.subplot(1, 3, 2)\n",
    "    plt.plot(history['temperature'], 'g-')\n",
    "    plt.xlabel('Iteration')\n",
    "    plt.ylabel('Temperature')\n",
    "    plt.title('温度下降曲线')\n",
    "    plt.grid(True, ls=\"--\")\n",
    "\n",
    "    # 子图3：k参数演化\n",
    "    plt.subplot(1, 3, 3)\n",
    "    plt.plot(history['k'])\n",
    "    plt.xlabel('Iteration')\n",
    "    plt.ylabel('k参数值')\n",
    "    plt.title('k参数变化曲线')\n",
    "    plt.legend(bbox_to_anchor=(1.05, 1), loc='upper left')\n",
    "    plt.grid(True, ls=\"--\")\n",
    "\n",
    "    plt.tight_layout()\n",
    "    plt.show()"
   ],
   "id": "ee66a50ae1eeb76d",
   "execution_count": 15,
   "outputs": []
  },
  {
   "metadata": {
    "ExecuteTime": {
     "end_time": "2025-08-27T13:46:40.253277800Z",
     "start_time": "2025-08-27T13:11:25.106835Z"
    }
   },
   "cell_type": "code",
   "source": [
    "def P(k_6):\n",
    "\n",
    "        # 定义一个方程组（微分方程组）\n",
    "    def pfun(y, t):\n",
    "        y0, y1, y2, y3 = y\n",
    "        return np.array([\n",
    "            1 / (m_A + delta_m) * (-k_3 * y0 - ro * g * s * y1 + f * np.cos(omiga * t) + k_5 * (y3 - y1) + k_6 * (y2 - y0)),\n",
    "            y0,\n",
    "            1 / m_B * (-k_5 * (y3 - y1) - k_6 * (y2 - y0) ),\n",
    "            y2\n",
    "        ])\n",
    "    t = np.arange(0, 500, 0.2)  # 创建自变量序列\n",
    "    soli = odeint(pfun, [0, 0, 0, 0], t)\n",
    "    # print(soli[:, 0])\n",
    "    # plt.figure(figsize=(5, 2))\n",
    "    # plt.plot(t, soli[:, 1], 'b', label=\"浮子位移\")\n",
    "    # plt.plot(t, soli[:, 3], 'r', label=\"振子位移\")\n",
    "    # plt.legend()\n",
    "    # plt.show()\n",
    "\n",
    "    re = 0\n",
    "    for i in range(200, len(soli[:, 2])):\n",
    "        re += k_6 * (soli[:, 0][i] - soli[:, 2][i])**2*0.2\n",
    "    return re/460\n",
    "re = P(40000)\n",
    "print(re)\n",
    "\n"
   ],
   "id": "a93479dfcaa0cc8b",
   "execution_count": 19,
   "outputs": []
  },
  {
   "metadata": {
    "ExecuteTime": {
     "end_time": "2025-08-27T13:46:40.255276900Z",
     "start_time": "2025-08-27T13:40:49.805539Z"
    }
   },
   "cell_type": "code",
   "source": [
    "from tqdm import tqdm\n",
    "\n",
    "y = []\n",
    "k = [i for i in range(1, 100100, 100)]\n",
    "for i in tqdm(range(1, 100100, 100)):\n",
    "    y.append(P(i))\n",
    "print(y)\n",
    "plt.figure(figsize=(4, 4))\n",
    "plt.title('平均输出功率与阻尼系数曲线')\n",
    "plt.plot(k, y, label='P(k)', )\n",
    "plt.xlabel('阻尼系数')\n",
    "plt.ylabel('平均输出功率')\n",
    "plt.legend()\n",
    "plt.grid(True, which=\"both\", ls=\"--\")\n",
    "plt.show()"
   ],
   "id": "8c9f89e7d3c0ed83",
   "execution_count": 46,
   "outputs": []
  },
  {
   "metadata": {
    "ExecuteTime": {
     "end_time": "2025-08-27T13:46:40.255276900Z",
     "start_time": "2025-08-27T09:57:01.208824Z"
    }
   },
   "cell_type": "code",
   "source": [
    "def SA(iter, t0, tf, alpha):\n",
    "    t = t0\n",
    "    kc = 20000\n",
    "    Pc = P(kc)\n",
    "    kb = kc\n",
    "    Pb = Pc\n",
    "    for i in range(iter):\n",
    "        kn = kc + np.random.randint(-50000, 50000)\n",
    "        kn = np.clip(kn, 0, 100000)\n",
    "        Pn = P(kn)\n",
    "        if Pn > Pc or np.random.rand() < np.exp((Pn-Pc)/t):\n",
    "            kc = kn\n",
    "            Pc = Pn\n",
    "            if Pc > Pb:\n",
    "                kb = kc\n",
    "                Pb = Pc\n",
    "        t = t*alpha\n",
    "        if t < tf:\n",
    "            break\n",
    "        print(f\"第{i}轮,kb:{kb},Pb:{Pb}, kc:{kc},Pc:{Pc}\")\n",
    "        history['temperature'].append(t)\n",
    "        history['best_P'].append(Pb)\n",
    "        history['current_P'].append(Pc)\n",
    "        history['k'].append(kb)\n",
    "    return kb, Pb\n",
    "re = SA(1000, 800, 0.001, 0.9)\n",
    "print(re)\n",
    "plot_SA_history()"
   ],
   "id": "5230104b2da26a12",
   "execution_count": 17,
   "outputs": []
  },
  {
   "metadata": {
    "ExecuteTime": {
     "end_time": "2025-08-27T13:46:40.283292900Z",
     "start_time": "2025-08-27T09:58:24.903919Z"
    }
   },
   "cell_type": "code",
   "source": [
    "for key in history.keys():\n",
    "    history[key].clear()\n",
    "print(history['temperature'])"
   ],
   "id": "c1670b45ca025376",
   "execution_count": 18,
   "outputs": []
  },
  {
   "metadata": {
    "ExecuteTime": {
     "end_time": "2025-08-27T13:46:40.286290200Z",
     "start_time": "2025-08-17T08:55:47.467594Z"
    }
   },
   "cell_type": "code",
   "source": [
    "# 在SA函数外部初始化记录容器\n",
    "history2 = {\n",
    "    'temperature': [],\n",
    "    'best_P': [],\n",
    "    'current_P': [],\n",
    "    'k': [],\n",
    "    'alpha': []\n",
    "}"
   ],
   "id": "99ec20b8d5701ec2",
   "execution_count": 13,
   "outputs": []
  },
  {
   "metadata": {
    "ExecuteTime": {
     "end_time": "2025-08-27T13:46:40.290290Z",
     "start_time": "2025-08-17T08:56:06.331558Z"
    }
   },
   "cell_type": "code",
   "source": [
    "def plot_SA_history2():\n",
    "    plt.figure(figsize=(10, 8))\n",
    "\n",
    "    # 子图1：损失函数变化\n",
    "    plt.subplot(2, 2, 1)\n",
    "    plt.plot(history2['best_P'], 'r-', label='Best P')\n",
    "    plt.plot(history2['current_P'], 'b--', alpha=0.5, label='Current P')\n",
    "    plt.xlabel('Iteration')\n",
    "    plt.ylabel('P')\n",
    "    plt.title('P变化曲线')\n",
    "    plt.legend()\n",
    "    plt.grid(True, which=\"both\", ls=\"--\")\n",
    "\n",
    "    # 子图2：温度衰减曲线\n",
    "    plt.subplot(2, 2, 2)\n",
    "    plt.plot(history2['temperature'], 'g-')\n",
    "    plt.xlabel('Iteration')\n",
    "    plt.ylabel('Temperature')\n",
    "    plt.title('温度下降曲线')\n",
    "    plt.grid(True, ls=\"--\")\n",
    "\n",
    "    # 子图3：k参数演化\n",
    "    plt.subplot(2, 2, 3)\n",
    "    plt.plot(history2['k'])\n",
    "    plt.xlabel('Iteration')\n",
    "    plt.ylabel('k参数值')\n",
    "    plt.title('k参数变化曲线')\n",
    "    plt.legend(bbox_to_anchor=(1.05, 1))\n",
    "    plt.grid(True, ls=\"--\")\n",
    "\n",
    "        # 子图4：alpha参数演化\n",
    "    plt.subplot(2, 2, 4)\n",
    "    plt.plot(history2['alpha'], c='orange')\n",
    "    plt.xlabel('Iteration')\n",
    "    plt.ylabel('alpha参数值')\n",
    "    plt.title('alpha参数变化曲线')\n",
    "    plt.legend(bbox_to_anchor=(1.05, 1))\n",
    "    plt.grid(True, ls=\"--\")\n",
    "\n",
    "    plt.tight_layout()\n",
    "    plt.show()"
   ],
   "id": "398b619d4091fcaf",
   "execution_count": 15,
   "outputs": []
  },
  {
   "metadata": {
    "ExecuteTime": {
     "end_time": "2025-08-27T13:46:40.293805Z",
     "start_time": "2025-08-27T13:31:25.357769Z"
    }
   },
   "cell_type": "code",
   "source": [
    "def P2(k, alpha):\n",
    "        # 定义一个方程组（微分方程组）\n",
    "    def pfun(y, t):\n",
    "        y0, y1, y2, y3 = y\n",
    "        return np.array([\n",
    "            1 / (m_A + delta_m) * (-k_3 * y0 - ro * g * s * y1 + f * np.cos(omiga * t) + k_5 * (y3 - y1) + k*(abs(y2-y0))**alpha * (y2 - y0)),\n",
    "            y0,\n",
    "            1 / m_B * (-k_5 * (y3 - y1) - k*(abs(y2-y0))**alpha * (y2 - y0) ),\n",
    "            y2\n",
    "        ])\n",
    "    t = np.arange(0, 500, 0.2)  # 创建自变量序列\n",
    "    soli = odeint(pfun, [0, 0, 0, 0], t)\n",
    "    # print(soli[:, 0])\n",
    "    # plt.figure(figsize=(5, 2))\n",
    "    # plt.plot(t, soli[:, 1], 'b', label=\"浮子位移\")\n",
    "    # plt.plot(t, soli[:, 3], 'r', label=\"振子位移\")\n",
    "    # plt.legend()\n",
    "    # plt.show()\n",
    "\n",
    "    re = 0\n",
    "    for i in range(200, len(soli[:, 2])):\n",
    "        re += k*(abs(soli[:, 0][i] - soli[:, 2][i]))**alpha * (soli[:, 0][i] - soli[:, 2][i])**2*0.2\n",
    "    return re/460\n",
    "re = P2(40000, 1)\n",
    "print(re)"
   ],
   "id": "a6ae873eb41a9ad",
   "execution_count": 36,
   "outputs": []
  },
  {
   "metadata": {
    "ExecuteTime": {
     "end_time": "2025-08-27T13:46:40.295804800Z",
     "start_time": "2025-08-27T13:38:25.435733Z"
    }
   },
   "cell_type": "code",
   "source": [
    "\n",
    "\n",
    "\n",
    "# 创建数据点 - 增加分辨率以获得更光滑的表面\n",
    "alpha_values = np.linspace(0, 100000, 10)  # 从0到100000，500个点（更密集）\n",
    "beta_values = np.linspace(0, 1, 10)        # 从0到1，200个点\n",
    "\n",
    "# 生成网格坐标\n",
    "Alpha, Beta = np.meshgrid(alpha_values, beta_values, indexing='ij')\n",
    "\n",
    "#计算Z值\n",
    "Z = []\n",
    "for i in tqdm(alpha_values):\n",
    "    for j in beta_values:\n",
    "        Z.append(P2(i, j))\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "from mpl_toolkits.mplot3d import Axes3D\n",
    "\n",
    "X, Y = np.meshgrid(alpha_values, beta_values) # X, Y 现在是 4x4 的二维矩阵\n",
    "Z = np.array(Z).reshape(X.shape)\n",
    "\n",
    "# 如果您自己的 z 数据是一个一维列表，需要根据网格大小重新塑形\n",
    "# 假设 z_list 是您的数据，并且顺序与网格点 (X.ravel(), Y.ravel()) 一一对应\n",
    "# Z = np.array(z_list).reshape(X.shape) # 这是最关键的一步！\n",
    "\n",
    "# 开始绘图\n",
    "fig = plt.figure(figsize=(10, 8))\n",
    "ax = fig.add_subplot(111, projection='3d')\n",
    "\n",
    "# 绘制曲面图\n",
    "# rstride 和 cstride 表示行和列的采样步长，越小越精细\n",
    "# cmap 是颜色映射方案，可选 'viridis', 'plasma', 'coolwarm', 'jet' 等\n",
    "surf = ax.plot_surface(X, Y, Z, cmap='coolwarm', edgecolor='none', alpha=0.8)\n",
    "\n",
    "# 添加颜色条\n",
    "fig.colorbar(surf, ax=ax, shrink=0.5, aspect=10)\n",
    "\n",
    "# 设置坐标轴标签\n",
    "ax.set_xlabel('X Label')\n",
    "ax.set_ylabel('Y Label')\n",
    "ax.set_zlabel('Z Label')\n",
    "\n",
    "# 设置标题\n",
    "ax.set_title('3D Surface Plot')\n",
    "\n",
    "# 显示图形\n",
    "plt.show()\n"
   ],
   "id": "2c6c113509214b73",
   "execution_count": 45,
   "outputs": []
  },
  {
   "metadata": {
    "ExecuteTime": {
     "end_time": "2025-08-27T13:46:40.297808900Z",
     "start_time": "2025-08-21T08:07:53.829065Z"
    }
   },
   "cell_type": "code",
   "source": [
    "def SA(iter, t0, tf, alpha):\n",
    "    t = t0\n",
    "    kc = 50000\n",
    "    ac = 0.5\n",
    "    Pc = P2(kc, ac)\n",
    "    kb = kc\n",
    "    ab = ac\n",
    "    Pb = Pc\n",
    "    for i in range(iter):\n",
    "        kn = kc + np.random.randint(-50000, 50000)\n",
    "        kn = np.clip(kn, 0, 100000)\n",
    "        an = ac + np.random.normal(0,0.5)\n",
    "        an = np.clip(an, 0, 1)\n",
    "        Pn = P2(kn, an)\n",
    "        if Pn > Pc or np.random.rand() < np.exp((Pn-Pc)/t):\n",
    "            kc = kn\n",
    "            Pc = Pn\n",
    "            ac = an\n",
    "            if Pc > Pb:\n",
    "                ab = ac\n",
    "                kb = kc\n",
    "                Pb = Pc\n",
    "        t = t*alpha\n",
    "        if t < tf:\n",
    "            break\n",
    "        print(f\"第{i}轮,kb:{kb},ab:{ab}, Pb:{Pb}, kc:{kc},ac:{ac},Pc:{Pc}\")\n",
    "        history2['temperature'].append(t)\n",
    "        history2['best_P'].append(Pb)\n",
    "        history2['current_P'].append(Pc)\n",
    "        history2['k'].append(kc)\n",
    "        history2['alpha'].append(ac)\n",
    "    return kb, ab, Pb\n",
    "re = SA(1000, 2000, 0.001, 0.95)\n",
    "print(re)\n",
    "plot_SA_history2()"
   ],
   "id": "273ca2ca10762f6b",
   "execution_count": 11,
   "outputs": []
  },
  {
   "metadata": {
    "ExecuteTime": {
     "end_time": "2025-08-27T13:46:40.297808900Z",
     "start_time": "2025-08-17T04:00:44.594343Z"
    }
   },
   "cell_type": "code",
   "source": [
    "for key in history2.keys():\n",
    "    history2[key].clear()\n",
    "print(history2['temperature'])"
   ],
   "id": "e707762b7dd7ce94",
   "execution_count": 141,
   "outputs": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
